Drag force on an impurity below the superfluid critical velocity in a 
quasi-one-dimensional Bose-Einstein condensate 
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The existence of frictionless flow below a critical velocity for obstacles moving in a superfluid is 
well established in the context of the mean-field Gross-Pitaevskii theory. We calculate the next 
order correction due to quantum and thermal fluctuations and find a non-zero force acting on a 
delta- function impurity moving through a quasi-one-dimensional Bose-Einstein condensate at all 
subcritical velocities and at all temperatures. The force occurs due to an imbalance in the Doppler 
shifts of reflected quantum fluctuations from either side of the impurity. Our calculation is based on 
a consistent extension of Bogoliubov theory to second order in the interaction strength, and finds 
new analytical solutions to the Bogoliubov-de Gennes equations for a gray soliton. Our results raise 
questions regarding the quantum dynamics in the formation of persistent currents in superfluids. 



Quantum fluids and their macroscopic manifestations 
have long intrigued physicists. Recent rapid advances 
in the experimental manipulation of dilute ultra cold 
gases 0, O, IS 3 S have opened a new window to 
the understanding of macroscopic quantum phenomena. 
Such experiments provide an exciting opportunity to ad- 
vance our theoretical understanding of quantum fluids as 
they probe regimes where calculations are more tractable 
than many traditional condensed matter systems. 

One of the most remarkable features is the emergence 
of superfluidity and the resulting critical velocity, below 
which an obstacle is able to move without experiencing 
any friction 0]. The presence of a Bose-Einstein con- 
densate (BEC) is commonly believed to be intrinsically 
related to superfluidity. In this paper we probe the su- 
perfluidity of a quasi- ID BEC with a delta- function im- 
purity moving through it at constant velocity. We in- 
troduce a complete set of zero eigenvalue solutions to 
the Bogoliubov-de Gennes equations for the excitation 
spectrum of the BEC and impurity. Our calculation is 
consistent to second order in the interaction strength. 
We find a drag force is present at subcritical speeds for 
all temperatures. We emphasize that this result is not 
in conflict with Landau's argument, as the drag force in 
this situation arises from the scattering of fluctuations 
(either quantum or thermal) rather than the creation of 
quasiparticles [7] . This mechanism was first suggested as 
a source of dissipation at zero temperature via a pertur- 
bative calculation for a uniform three-dimensional BEC 
using the Born approximation [8[ . The current work goes 
beyond the perturbative regime and is applicable to re- 
pulsive delta-function impurities of any strength. Fur- 
thermore we calculate the force for systems at finite tem- 
perature, and find a clear distinction between the quan- 
tum and thermal regime. These results could be tested 
experimentally using foreseeable technology . 
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where $ = ^(x,t) is the bosonic field operator obey- 
ing the usual equal- time commutation relations, and 
g > is related to the 3D s-wave scattering length a: 
g ~ 2ftw±a (we have assumed a sufficiently tight trans- 
verse trapping potential of frequency co± [10]). We use 
the following dimensionless coordinates for the system. 
Time: r = (gnoo/h)t, and length: z = x/£ where 
£ = h/ y/mgrioc is the healing length, and is the 
total density, at x = ±oo (far away from the impu- 
rity). The impurity strength and impurity velocity are 
parameterised by fj = rj/(v s H) > and v — v/v s re- 
spectively, where v s = ygn^o/m is the speed of sound 
far from the impurity. Finally we rescale the field op- 
erator, j/j = ^/riooif. In these units the commutation 
relations become [<p(z, r), 0^(z', r)] = ^S(z — z') where 
7 = mg I \h 2 noo) is the ratio of interaction energy to ki- 
netic energy. Our calculations are based on the assump- 
tion that e = 7 1 / 4 <C 1 and this defines the small pa- 
rameter in our perturbative expansion. The equation of 
motion for (p is then 
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An understanding of the relevant length scales in a ID 
Bose gas is important [11[. The ground state of an in- 
finitely extended, weakly interacting ID Bose gas is that 
of a quasicondensate [12], with phase coherence length: 
= £ exp( v // 7r 2 /7). Hence phase coherence is present 
across regions of size 1$ ^> £, but not over all space. We 
define the size of our system, L, to be the same order as 
the phase coherence length. That is L ^> £ and L < Z^, 
such that we have a true condensate which can be con- 



The second-quantised Hamiltonian of the system in the sidered homogeneous [24( . We emphasize that we are not 



2 



considering regimes of ID systems for which there does 
not exist a true-BEC [13[ . 

We ultimately wish to calculate the force on the im- 
purity, given by the gradient of the potential, F = 
— {(p^d z [fjS(z)] 0). First we find the solution to Eq. ([2]) 
from Bogoliubov's perturbation expansion (p ~ ((po + 
e(fi +e 2 (/?2)e _ ^^ /i0+e/il+e The expansion is conceptu- 
ally summarised as (i) cpo is the condensate wave function 
in the absence of all fluctuations; (ii) (p\ describes the 
fluctuations of the field in the presence of a condensate 
cpo; (hi) cf2 is the modification of the condensate wave 
function due to the presence of fluctuations. 

(i) O(e ): At zeroth order the commutators vanish 
and a c-number field describes the system. The con- 
densate wave function is given by the solution to the 
Gross-Pitaevskii equation 



L[(fo](fo = 0, 



(3) 



where L[(p ] = -\d 2 +ivd z + fj6(z) + \(fo\ 2 - 1 (/i = !)• 
This can be solved analytically [14[ for impurity veloci- 
ties less than the critical velocity v < v Cl where v c de- 
pends on fj by v 2 c = 1 - 2 -f + ^ (| - R) + rf (f + ^) 

where R 3 = 2y/2f\/ (yi + 87f/27 - l) 0. To do so, 

the term fjS(z) in Eq. (|3|) is replaced with a derivative 
jump: d z cpo(z = + ) — d z (po(z = 0~) = 2f]. The dark 
soliton solutions [17] are used to give 



e ia ~ {cos(0) tanh [z~] 
e ia + {cos(0) tanh [z+] 



-ism(0)} (z < 0), 
-ism(0)} (z>0), 



(4) 



where the impurity velocity is parameterised by v = 
sin(0) < v C: and zf = (z ± zq) cos(0). The phases <7 + = 
— 0, and cr_ = — ir + 2 arctan[sin(2#) / (exp(2zo cos 6) — 
cos(20))] are fixed, although cpo has an arbitrary global 
phase factor. Here, we work in the limit of a large 
system, where the effects of fluctuations in the to- 
tal number of particles are negligible [l8[ • The quan- 
tity zq is determined by the derivative jump con- 
dition, which reduces to numerically solving fj = 
cos 3 (#) tanh [cos(6)z ] / (sin 2 ((9) + sinh 2 [cos(6)z ]) for a 
given value of fj and v < y a . Two solutions exist for 
zq but only one is physical [14]. Figure [T] shows a plot of 
<fo(z) f° r a specific fj and v. As expected, at this level 
of approximation the drag force F = 0. For numerical 
solutions to Eq. (|3]) in higher dimensions see Ref. [15[. 

(ii) 0(e 1 ): The quantum and thermal depletion of the 
condensate is accounted for at first order. Linearisation 
of Eq. ([2]) with respect to e results in 



id T (pi = L[V2(p ] 
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The density is unaffected at order e (as (0i) = 0) and 
therefore \i\ = 0. A Bogoliubov transformation is then 

performed, (p\ = jdk Uk{z)e~ %EkT c\k + v^{z)e %EkT c\\ 
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FIG. 1: (color online). Mean field solutions for (a) the con- 
densate amplitude |<£o(<z)| and (b) the condensate phase from 
Eq. J4j for fj — 1 and v/v s — 0.35. The dashed line shows the 
continuation of each of the underlying dark soliton solutions. 



where (a\) are the annihilation (creation) operators 
of the excitations (note that k is a continuous index). 
They obey the usual commutation relations, [o^,o^,] = 
5(k—k') and [a^a^/] = = 0. The quantum state 

of the system is defined such that the occupation number 
of each excitation is rik = {a\ak) = l/(e Ek ^ T — 1). The 
amplitudes, Uk and Vk, are determined by 



L[y/2<Po] 



9 

ft 



Uk — EkUk, 



(6) 



where Uk = (uk,Vk) T . The normalisation which pre- 
serves the bosonic commutation relations is given by 
f dzul(j z Uk> = 5(k — k') where o z is a Pauli matrix. 
The solutions of Eq. (j6]) for z G M ± are [H, 



(I 



k 



i cos 6 tanh (zf)Y 



(§-^ + icos0taiih(*±))' 



, (7) 



where Ek has two distinct branches: 
±k (rsin(<9) + 



Ek = 

1^ , corresponding to the left 
and right propagating excitations and hence Ek > 
(special attention must be paid to the Ek = mode). 
For a fixed excitation energy E, there are four possi- 
ble wavenumbers, given by the four roots of (l/4)/c 4 + 
cos 2 (6)k 2 + 2Es'm(6)k - E 2 = 0. Two of the roots will 
be in R, denoted k and k r with one positive and one 
negative. k r is a reflected mode, and k r = —k only 
when v = 0. For u 7^ a Doppler shift will occur. 
The other two roots will be complex conjugates, (with 
nonzero imaginary parts) denoted k e and fc*. On either 
side of the impurity, one must select the exponentially de- 
creasing mode to satisfy the boundary conditions of the 
problem. Each wavenumber corresponds to a particu- 
lar solution of Eq. (|6]), which is ultimately a 4th order 
differential equation. The normalisation of the eigen- 
vectors in Eq. is N 2 J dzQa z ^k = S(k - k') where 



\^kJ^(k 2 /A + l) 



-1/4^-1 



Thus the solutions in 



Eq. (O correspond exactly to the usual Bogoliubov ex- 
citations of a uniform BEC traveling at velocity v for 

z — > ±00 



21|. 



To incorporate the impurity at z = 0, we separate the 
problem into two independent scattering problems: One 
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in which the incoming mode approaches from z = — oo 
and the other in which it approaches from z = +00. The 
boundary conditions we impose are: (i) the amplitude of 
the incoming wave is given by N k in accordance with the 
assumption that, away from the impurity, the BEC has 
no knowledge of the impurity. Implicit in this condition 
is a time scale over which it is assumed the scattered 
waves have not reached edge of the system, (ii) Scattered 
waves cannot exponentially increase as they move away 
from the impurity, (iii) Scattered waves must be causally 
related to the incoming wave. For the case where the 
incoming wave approaches from z = — 00 the solution is 




rtk* (*<0), 

(z > 0), 



(8) 



where k > and the complex constants A r , A tl A etl and 
A er are completely determined by u k \ z = + = u k \ z = - 
and d z Uk\ z =o+ — d z u k \ z = - = %Vu k \ z =o, which arise from 
the fjS(z) term in Eq. (|6]). Solving the scattering prob- 
lem where the incoming wave approaches from z = +00 
follows the same logic, except now k < and the z > 
and z < conditions in Eq. JHJ) have to be swapped. 

At this stage the problem has been reduced to a set 
of linear algebraic equations for each value of k. Ana- 
lytic progression is hindered by the fourth order polyno- 
mial determining k e and k r . However, quantities such 
as: (<p{(z)<pi(z)) = J dk[\u k \ 2 n k + \v k \ 2 (l + n k )} and 
((fi(z)(fi(z)) = J dk[u k vl(l + 2n k )] are easily attained 
numerically to very high accuracy once one fixes 77, v, and 
T. One must include a small k cut-off [of order ~ 1/L] to 
prevent long wavelength fluctuations destroying the long 
range order in the system. Results are logarithmically 
dependent on the choice of this cut-off [22I ]. 

(iii) 0(e 2 ): At second order, the generalised-Gross- 
Pitaevskii equation [l8[ is 



H\ V2 {z)) = \f(z)), 



(9) 



where f(z) = 2cp J dk[\u k \ 2 n k + |u fe | 2 (n fc + 1)] + 
<Pq J dk[ukV^,(l + 2rife)] — [12^0, (this definite integral is 
performed numerically using the analytic results of the 
previous section), and 



n 
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(10) 



Here, a ket always corresponds to a complex conjugate 
pair: |#) = (• , •*) T . The shift in the chemical poten- 
tial, 112 is defined on either side of the impurity such 
that f(z) — > as z — > ±00, and contains the necessary 
terms to ensure orthogonality between the condensate 
mode and the excitations [181 ]. 

To solve Eq. |9|) for ^2(2), we use four linearly indepen- 
dent solutions of the homogeneous equation, H\uo) = |0), 
to construct a 2x2 Green's matrix, Q{z, s). Finding these 
four linearly independent solutions constitutes a nontriv- 
ial step, and we were unable to find expressions for them 



in the previous literature; we provide a detailed presenta- 
tion of these solutions in [23[ . Briefly, linear combinations 
of Eqs. (2-5) of [23[ are used to generate four solutions, 
denoted |Aj(s)), appropriate to the presence of an im- 
purity [the f]5(z) term in Eq. (|SJ)]. The Green's matrix 
satisfies 



Hg(z,s)=6(z-s)I 2 , 



(11) 



where I2 is the 2x2 identity. This condition yields 
Eqs. (13-16) of Ref. [23!]. The boundary conditions for 
Eq. ([9j) are: (i) cp 2 does not exponentially increase away 
from the impurity; (ii) cp 2 is symmetric about the im- 
purity for v = 0; (iii) Any non-zero drag force on the 
impurity acts to decrease the relative motion between 
condensate and impurity. These conditions combined 
with the symmetry of the Green's matrix [specifically 
Q(z, s) = G^(s, z) which arises from the operator TL being 
Hermitian] uniquely define the solution 



M*))= dsg(z,s)\f(s)), 



(12) 



up to a global phase factor. We note that Eq. (fT2|) is 
not a general formula for any function f(s) — the fact 
that f(s) decays to zero as s — > ±00 compensates for the 
linear divergences in Q. The integral in Eq. ([T2]) is also 
performed numerically. 

We now evaluate the force to 0(e 2 ) 



F~fjd z 



dk[\u k \ 2 n k + \v k \ 2 (l +n k )] +2Re{(/? ^2} 



(13) 

where all the derivatives are performed analytically us- 
ing Eq. dSj) and Eqs. (2-5) of [23[. The results are shown 
in Fig. [Sfa) for the case of zero temperature. The fact 
that the force is nonzero below the critical velocity may 
be unexpected, although we note the magnitude is much 
smaller than the supercritical forces present above v c [3] • 
We emphasise that quasiparticle creation is not the mech- 
anism behind this drag force, and consequentially this is 
not a contradiction of Landau's statements 0. At zero 
temperature fluctuations in the quantum vacuum scat- 
ter from each side of the impurity, and the asymmetric 
Doppler shift in the reflected waves results in an imbal- 
ance and a non-zero net force. A crucial assumption in 
this result is that the incoming waves have no knowledge 
of the impurity. In any finite system this assumption 
would only be valid up to a time t c ~ L/v s after the 
initial acceleration of the impurity. It could be imagined 
that the force on the impurity will decay to zero over this 
time scale t c , which can be thought of as the relaxation 
time of the quantum vacuum. We intend to test this 
conjecture by performing dynamical simulations that in- 
corporate the effects of the quantum fluctuations; this 
significant project lies beyond the scope of this work. 

The results at finite temperature are shown in Fig.[2|(b- 
d) . The quantum and thermal regimes are clearly distin- 
guished with a cross-over near T ~ 10 -3 gn^. The force 




10" 5 10" 4 10" 3 10" 2 10" 1 10" 3 

T/gn^ T/gn 



FIG. 2: (color online). The drag force on the impurity (in 
units y/jgn^). (a) Drag force at T = as a function of 
velocity, fj — 0.3, 1,3 for diamond, square, circle respectively, 
and each line terminates at their respective critical velocities 
Vc/vs = 0.64,0.37,0.16. (b-d) Drag force as a function of 
temperature, (b) 77 = 0.3 and velocities: blue diamond v = 
0.14, blue square v — 0.35, blue circle v — 0.56. (c) fj = 1 
and velocities: red diamond v — 0.08, red square: v — 0.2, 
and red circle v = 0.3. (d) is data taken from (b) and (c), 
and compares the drag on a strong impurity moving at slow 
speed with that on a weak impurity moving at high speed. 
The strong impurity can have a larger drag in the quantum 
regime while the weak impurity will have the larger drag in 
the thermal regime. 

increases linearly with T in the thermal regime due to 
the scattering being most prominent for the low energy 
modes whose occupation is n& = l/(e Ek ^ T — 1) ~ T '/ ' 
for T ^> Ek- The reflected modes for the high energy 
incoming waves have a negligible Doppler shift and only 
a small contribution to the force. 

In conclusion, we have calculated the drag force acting 
on an impurity moving through a true BEC, in a quasi- 
1D geometry, using analytic solutions of the Bogoliubov 
perturbation expansion. We find that a nonzero force 
arises for zero and finite temperature at all velocities due 
to an imbalance in the Doppler shift in the scattering of 
quantum and thermal fluctuations. The crossover from 
the quantum to the thermal regime occurs at tempera- 
tures near T ~ 10~ 3 gn oo . We have proposed a mecha- 
nism that predicts a time scale of Lj v s for which it would 
be possible to observe the force in a finite system. 
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AUXILLIARY PUBLICATION FOR: DRAG FORCE ON AN IMPURITY BELOW THE SUPERFLUID 
CRITICAL VELOCITY IN A QUASI-ONE-DIMENSIONAL BOSE-EINSTEIN CONDENSATE 

Here we describe the solution of Eq. (11) in the main text (repeated here for convenience) 

ng(z,s) = 5(z-8)i 2 , (14) 

where I2 is the 2x2 identity. We begin by writing down the four linear independent solutions to the homogenous 
equation 

H\lo) = |0> (15) 

for the separate cases of z G M + and z G R~. These are denoted \uf(z)) (recall that the ket notation refers to a 
complex conjugate pair) where, 

uf(z G M ± ) = ipf, (16) 

uf(z G M ± ) = e ia± sech 2 (^), (17) 

uf(z G M ± ) = e ia± (sech 2 (z±) [2zf - zf cosh(2z c ± ) + (3/2) sinh(2^)] tan((9) + 2i [zf tanh(z c ± ) - l]) , (18) 
uf(z G M ± ) = e iCT± sech 2 (^) [zf (10 - 4cos 2 ((9) - 8sin(0) sin((9 - 2izf)) + 

cosh(z c ± ) [i sin (2(9 - 3izf) - U sin (2(9 - + 6 sinh(2^)] . (19) 

To the best of our knowledge, Eqs. ([15]) and ([T9]) do not appear in the previous literature. Note from these solutions, 
the behaviour as z — > ±00: — > constant, cj^ — > 0, cj^ — > linear increase/decrease, and uf — > exponential 
increase/decrease, which can be seen from Fig. [3] 

Using these eight functions \u^(z)) (where j always runs over the indices j = 1, ... ,4) we construct four linearly 
independent solutions (denoted \Aj(z))) to Eq. ([15]) for the entire domain z G R. This is done by matching the 
solutions in R ± across the point z = subject to the conditions 

\A j (z = + ))=\A j (z = 0-)), (20a) 
a,|A 5 -(z = 0+)) - d z \Aj(z =0")) = 2t7|A,(z = 0)), (20b) 




FIG. 3: (color online). The zero eigenvalue solutions of Eq. JT5J) for zq = and cr± = as a function of z and 6. (a) and (b) 
show the real and imaginary parts of uji(z), (c) and (d) show the real and imaginary parts of U2(z), (e) and (f) show the real 
and imaginary parts of u>3(z), and (g) and (h) show the real and imaginary parts of uj^z). 
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which are the usual conditions appropriate for the term rjS(z) in the differential operator TL The result is 

' Al W> = {l5w> (21) 

|A2W) = { + a 2 \cu+(z)) + a 3 \co+(z)) + a 4 \cu+(z)) (22) 

lAs{z)) = { £|j£(l)> + & 2 |u> 2 + (*)> + &3|<4(*)> + bM(z)) (23) 

|A4(Z)) = { £|Jf(l)> + c 2 |u> 2 + (*)> + c 3 | W +(z)) + c 4 K + (*)> (24) 



Kl(s)) 


- |AiW> 


K 2 (s)) 


- Ms)) 


K 3 (s)) 


- \Hs)) 


K 4 (s)) 


- |A 4 ( S )) 



(27) 
(28) 



where a^, bj,Cj G M and are determined from Eqs. (|2Q|) above. 

The next step involves finding Q(z, s) which is a solution to Eq. (11) of the main article. We write 

f |A 1 (^))(A 1 ( S )| + |A 2 (z))(A 2 ( s )| + |A 3 (z))(A300| + |A 4 (z))(A 4 ( s )| (z < s) 
S(z,s)={ (25) 
[ IAxWXkxWI + \A 2 (z))( K2 (s)\ + \A 3 (z))(k 3 (s)\ + \A 4 (z)) (k 4 (s)\ (z > s) 

With some work it can be shown that the conditions 

Q(z = s + , s) =G(z = s~ , s) (26a) 

= s+, s) - 0*0(2 = a", s) = -2I 2 (26b) 

are completely equivalent to the following equations 

isec 2 (0)|A 3 ( S )>, 

isec(0)tan(0)|A 3 ( S )> + ^ sec 3 (0)|A 4 ( S )), 

- X - sec 2 (0)|Ai(s)> - J sec(0) tan(0)|A 2 (s)), (29) 

__L sec 3 W |A 2(s)) . (30) 

This is a satisfying result that can be anticipated from the fact that the operator H is self- adjoint, and that vectors 
\Xj(s)) and \Kj(s)) should always be solutions to the adjoint problem. 
To remove the components of Q which have an exponential divergence at z = ±00, we impose the conditions 

Ms)) = 10), (31) 
a 4 \K 2 (s)) + b 4 \K 3 (s)) + c 4 \k 4 (s)) = |0). (32) 

Enforcing the symmetry of the Green's matrix Q(z, s) = Q^(s, z) results in 

4 = 0, (33) 

\t = 0, (34) 

nl = \l (35) 

k\ + ^-kI = - j sec(0) tan(0) (36) 

"A? = ±4, (37) 

where \Xj(s)) = £* =1 A«|A g ( S )) and \ Kj (s)) = k?|A,(s)>. 

The final condition to ensure the Green's matrix is unique relates to how the system responds to the broken 
symmetry caused by the moving impurity. The force could be either positive or negative, and the physically correct 
solution is when the force is positive as this is the only energy conserving solution. This gives 

A} = 0, (38) 
A 3 = 0, (39) 
A 2 = 0, (40) 
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FIG. 4: (color online). Shows the Green's matrix for 77 = 1 and v — 0.3. (a) and (b) show the real and imaginary components 
respectively of Gu(z, s). (c) and (d) show the real and imaginary components respectively of Gi2(z, s). 



and completes the solution to Eq. (11) of the main article. 



